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1.  INTRODUCTION 


In  any  estimation  process  there  are  four  quantities  which  may  have  to  be  taken  account  of : 

(i)  bias  of  estimates ; 

(ii)  variance  of  estimates; 

(iii)  time  taken  to  obtain  estimates  on  machine  that  is  available; 

(iv)  memory  and  storage  requirements  of  method. 

The  particular  problem,  for  which  the  following  method  was  generated,  is  to  estimate  the 
poles  of  a  linear  dynamical  system  in  close  to  real  time.  The  linear  system  is  the  representation 
of  an  aircraft  in  flight  and  the  poles  are  determined  for  stability  reasons. 

It  is  required  to  know  how  these  poles  vary  as  a  function  of  forward  speed  and  altitude. 

In  this  instance  the  aircraft  will  fly  at  a  constant  speed  for  a  period  of  time.  Five  minutes 
will  be  available  for  calculations  then  the  speed  will  be  incremented  or  altitude  changed. 

All  calculations  have  to  be  done  on  a  minicomputer  which  has  a  storage  capacity  of  16K 
and  can  calculate  a  2048-point  Fast  Fourier  Transform  (FFT)  in  approximately  1  -5  minutes. 

Any  non-linear  method  is  effectively  ruled  out  due  to  this  time  constraint. 

In  a  more  general  setting,  where  “time”  is  not  so  crucial  a  factor,  even  a  non-linear  method 
like  maximum  likelihood  or  least  squares  needs  some  initial  estimates.  This  method  will  give  good 
initial  estimates  for  any  non-linear  method. 

Although  the  title  of  this  note  may  imply  that  the  estimates  are  both  quick  and  inaccurate, 
this  is  not  the  case. 

Resolution  is  much  better  than  that  of  power  spectral  methods  as  are  estimates  of  dampings 
and  frequencies. 

Providing  that  the  maximum  number  of  frequencies  of  decay  which  are  “close”  together 
is  two,  then  the  greatest  level  of  mathematical  complexity  is  the  solving  of  two  simultaneous 
equations  and  the  solving  of  a  quadratic  equation  with  complex  coefficients. 

More  complicated  situations  require  the  solution  of  larger  sets  of  simultaneous  equations 
and  higher  order  polynomial  equations ;  the  essential  simplicity  of  the  method  is  then  lost  and 
it  is  then  preferable  to  go  to  a  full  maximum  likelihood  method.1 


2.  THE  PROBLEM 

Given  a  signal  ^,Aie~a*t  cos  (toft-tfa)  it  is  desired  to  estimate  wi  and  «i  where  cm  is  the  decay 
frequency  and  <xi  is  the  decay  rate. 

It  is  common  to  estimate  /(  and  ««/»<  where  /c  is  the  frequency  in  hertz  (2vft  =  «m)  and 
«(/<m  is  normally  expressed  as  a  percentage.  The  reason  for  this  choice  of  parameters  is  simply 
that  for  stability  reasons  we  would  like  the  poles  of  the  system  to  lie  within  a  certain  sector  in 
the  complex  domain.  Care  must  be  taken  with  this  concept  especially  for  low  frequencies. 


2.1  Single  Pole-Single  Zero  Approximation 
In  a  small  frequency  band  the  Fourier  transform  is  assumed  to  behave  as 

No+Ng 

Aj+r 


(2.1.1) 


where  *  =  ia>  and  No,  Ni,  D»  may  be  complex.  Conjugate  poles  are  ignored  as  neighbouring 
frequencies  will  have,  in  general,  a  greater  influence. 

Only  three  frequency  points  are  needed  to  obtain  No,  Nt  and  Da-  Guidelines  will  be  presented 
later  as  to  how  these  points  may  be  selected. 
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Let  the  three  frequency  points  be  si,  si,  ss  and  let  the  corresponding  values  of  the  transform 
be  Zi,  Zt,  Zt.  Equating  observed  values  with  expected  values  it  follows  that 

(iVo+A'isi)  =  Zi(s\+Do),  (2.1.2) 

(No+NiSi)  =  Z^st+  Do),  (2. 1 .3) 

(No+Nisa)  =  Zs(s3+  Do).  (2. 1 .4) 

Subtracting  (2.1.3)  from  (2.1.2)  and  dividing  by  (si— f2)  it  follows  that 

Ni=s^1~^  +  DgZ}~Z3 

S1—S2  Si— Si 

Performing  a  similar  operation  on  (2.1.3)  and  2.1.4)  it  follows  that 

Ni  =  SiZi-*Zi+DgZi-Zi  (2.1.6) 

Ss—  Ss  Si— si 


siZi  —stZi 


StZt  —ssZs 


Zi-Zt 


Zt-Zt 


Y\  and  Yt  are  the  so-called  divided  differences  of  the  function  sZ;  X,  and  Xt  are  the  divided 
differences  of  the  function  Z. 

Subtracting  (2.1.5)  from  (2.1.6)  and  solving  for  Do  it  follows  that 

Do  =  -  (  Yi  -  Yt)/(Xi  -  Xt).  (2. 1 .7) 

D0  contains  the  parameters  of  interest  and  is  of  the  form  «— iw.  If  required  Ni  and  No  may  be 
found  by  hack  substitution  in  equations  (2.1.6)  and  (2.1.4). 

2.2  Equation  of  Nyqoist  Locus 

The  complex  response  Z(iw)  may  be  rewritten  in  the  form 


Z(iw)  =  Ni  + 


Oo+  iV 


where  N  —  No—NiDo. 


Now  let  N  =  A+iB,  Do  =  a—iw,  Ni  =  Xo+iYo,  X(iw)+iY(im)  =  Z(iw).  Then 

v,.-  ^  V  ,  l) 

X(iw)  =  Xo+  J+^—yt- 


Y(iw)  =  Yo  + 


fla  a>i) 

«*  -f  (<o  — 


After  some  algebra  it  can  be  shown  that 


Equation  (2.2.3)  is  recognised  as  the  equation  of  the  circle:  centre  ^Xo  +  To  + 

\Po+B* 

radius  — - - -.  It  should  be  noticed  that  it  is  possible  to  have  a  poor  estimate  of  the  circle, 

2a 

which  has  parameters  involving  the  zero  as  well  as  the  pole,  whilst  still  obtaining  good  estimates 
of  the  pole.  That  is,  it  is  more  important  to  determine  how  ‘fast’  the  Fourier  transform  moves 
along  the  circle  than  to  define  the  circle  precisely. 
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13  Two  Zcros-Two  Poles  Approximation 

If  two  damped  frequencies  are  “close”  together  then  a  two  zero-two  pole  approximation 
is  called  for: 

-7/  \  No+Nis+Nis* 

2(1  =  T+lhs+Dv*'  (231) 

where  r  =  im.  Equation  (2.3.1)  has  five  complex  parameters  hence  five  frequency  points  are 
needed. 

Equating  expected  values  with  observed  values  gives  five  equations 

Zi{l +  !>,*+ />»*<*}  =  Afo+Afo+JV**2  i  =  1.2 . 5.  (2.3.2) 


„ , _  1  (Zt—Zu i  _  Zi+i—Zt+t) 

St  — St+i  S(+i  — S(+2j 

Y  (i)  —  1  (StZt— St+iZt+i  _  St+iZt+i—St+zZt+j) 

St— S<+2(  St—Si+i  Si+i—Si+i  j 

2 /j-\ _ 1  (SfZj—Sf+iZj+i  ~  Sj+iZi+i—Si+tZi+t) 

St— Sm}  Si-St  Si+i—Si+t  J 


Si— Sr+2^  Si-St  Si+i—Si+t  ) 

•*s(i),  K2O).  Z-^i)  are  the  second  divided  differences  of  Z,  sZ  and  s2Z  respectively. 
The  equations  (2.3.2)  can  now  be  reduced  to 


X^i)+DxY2(i)+DtZHi)  =  Ni 


If  we  now  let 


Xt  =  Xt(i)-X2(i+ 1);  Y,'  =  YHi)-  Y2(i+ 1);  Zt  =  Zs(i)-Z2(i+ 1)  *  =  1.2 

Then  (2.1.3)  reduces  to 

X,  +  Di  Y,'  +  D2Zt  =0  i'=l,2. 


„  XiYi-XiYY 
Dt~  Zi'  Yt'-ZtYY 

(2-3.5) 

XiZY-Xt'Zi 
1  -YiZ2'+  YYZY  , 

The  required  poles  are  now  the  solution  of  1  4-  Dis 4-  Dts1  =  0,  i.e. 

s  =  -  Oi±Vf)i2-4Os  (2.3.6) 

2  Os 

Di  *  0. 

It  should  be  noticed  that  (2.3.6)  involves  a  complex  square  root. 

No  attempt  will  be  made  to  determine  the  locus  in  the  Nyquist  Plot  as  the  parametric 
representation  (2.3.1)  is  probably  the  simplest.  Furthermore  it  is  not  the  curve  itself  but  how 
"fast”  the  Fourier  Transform  moves  along  the  curve  which  is  important.  Back  substitution  will 
determine  No,  N\  and  Nt  if  required. 


3.  SELECTION  OF  FREQUENCY  POINTS— SOME  GUIDELINES 

The  representations  of  the  response  are  only  local  representations  and  hence  do  not  cover  a 
wide  frequency  band.  Ideally  the  points  selected  should  cover  the  damped  frequency.  Further¬ 
more  as  successive  differences  are  calculated  only  points  with  large  signal-to-noise  ratios  should 
be  used.  From  equation  (2.1.7)  (single-pole  case) 

Do  —  -(Yi-Y^UXi-Xn). 
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To  reduce  sensitivity  it  is  necessary  that  X1—X2  is  large,  i.e. 

Zi-Zn  ^  Zj-Zj 
S1—S2  st—sa 

is  large.  For  equally-spaced  points  this  means  searching  for  local  maxima  of  |Zi— 2Zj+Zj|. 

Interpolated  points  in  a  Fourier  transform  should  not  be  used  as  these  will  not  be  statistically 
independent. 

Interpolated  points  in  a  Fourier  transform  come  from  using  the  shifting  theorem,  i.e.  by 
taking  the  FFT  of  [A"(f  V-*"*]  (a>  small)  rather  than  the  FFT  of  X(t). 


4.  DISCUSSION 

The  method  presented  was  devised  because  traditional  methods  (e.g.  half-power  points  and 
rapid  phase  changes)  have  a  number  of  disadvantages. 

The  disadvantages  of  the  half-power  method  are  that: 

(i)  it  uses  only  information  from  the  power  spectrum  and  ignores  the  phase  spectrum; 

(ii)  the  peak  power  may  in  fact  be  missing  as  may  the  half-power  points  since  the  FFT  is 
discrete; 

(iii)  the  method  does  not  allow  for  the  presence  of  a  nearby  zero  in  which  case  the  peak 
power  does  not  occur  at  the  frequency  of  decay. 

A  nearby  zero  can  occur  quite  frequently  especially  if  the  transfer  function  is  minimum  phase 
(the  response  is  measured  quite  close  to  the  excitation  point). 

Similar  criticisms  may  be  levelled  against  methods  which  use  only  phase  information. 

The  method  devised  in  this  note,  which  uses  both  amplitude  and  phase  information,  allows 
up  to  two  zeros  and  two  poles  to  be  close  together.  The  stability  of  the  estimates  may  be  checked 
by  seeing  if  points  other  than  those  used  in  the  estimation  lie  near  the  Nyquist  plot. 

The  stability  of  the  estimates  may  also  be  checked  by  using  a  different  selection  of  points. 

The  situation  may  arise  where  it  is  not  dear  if  there  are  any  poles  in  a  certain  frequency  band 
or  not.  The  stability  of  the  estimate  is  then  quite  important  for  deciding  an  answer  to  this  question. 

It  has  been  assumed  throughout  that  the  frequencies  of  interest  are  well  away  from  the 
Nyquist  frequency.  By  using  the  Fourier  sum  rather  than  the  Fourier  integral  it  is  possible  to 
change  the  method  so  as  to  work  near  the  Nyquist  frequency. 


5.  EXAMPLES 

Figures  la,  2a,  3a  are  the  time  history,  amplitude  spectrum  and  Nyquist  plot  respectively  of 
a  signal  without  noise.  The  time  history  is  generated  mathematically  and  contains  decay  frequen¬ 
cies  of  20  Hz,  55  Hz  and  60  Hz  with  decay  rates  of  2,  1  and  5%  of  the  corresponding  decay 
frequency  respectively. 

A  single-pole  curve  is  fitted  to  points  at  17-5  Hz,  20  Hz  and  22-5  Hz  (the  position  of  the 
first  point  is  indicated  by  a  “1”  in  Figure  3a).  The  decay  frequency  and  damping  as  determined 
by  equation  (2.1.7)  are  20  Hz  and  2-02%. 

A  two-pole  curve  is  fitted  to  points  at  52-5,  55,  57-5,  60  and  62-5  Hz  (the  position  of  the 
first  point  is  indicated  by  a  "2”  in  Figure  3a).  The  two  decay  frequencies  as  determined  from 
equation  (2.3.6)  are  55  and  60-01  Hz,  and  the  dampings  are  1  and  5-04%  respectively. 

Figures  16,  2b,  36  are  the  corresponding  plots  with  noise  added  to  the  signal.  The  single-  and 
double-pole  curves  are  fitted  to  the  same  frequency  points  as  for  the  original  signal. 

The  r.m.s.  value  of  the  noise  is  Aly/l  where  A  =  J  largest  |Amplitude[.  This  implies  that  the 
signal- to-noise  ratio  is  always  <5>/3.  The  estimated  decay  frequencies  are  20,  55-09  and  61  -4 
Hz  while  the  corresponding  dampings  are  1-8,  1-11  and  5-07%. 
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f  6.  CONCLUSION 

~  simple  procedure  for  determining  initial  estimates  of  decay  frequencies  from  Nyquist 
plots  has  been  determined.  The  process  gives  acceptable  results  even  in  the  presence  of  a  moderate 
amount  of  noise.  These  initial  estimates  may  then  be  refined  by  using  some  non-linear  method 
if  required^ 
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FIG.  3(b)  NYQUIST  PLOT  FROM  TIME  HISTORY  WITH  NOISE 
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